\documentclass{article}
% Packages --------------------------------------------------------
\usepackage[margin = 0.8in]{geometry}
\usepackage{amsmath,amssymb,amsthm,amsfonts}
\usepackage{graphicx}
%\usepackage{hyperref}
% Theorems -------------------------------------------------------
\newtheorem{thm}{Theorem}[section]
\newtheorem{cor}[thm]{Corollary}
\newtheorem{lem}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}
\theoremstyle{definition}
\newtheorem{defn}[thm]{Definition}
\theoremstyle{remark}
\newtheorem{remark}{Remark}
\newtheorem{rem}[thm]{Remark}
\numberwithin{equation}{section}
\newcommand{\var}{\operatorname{var}}
\newcommand{\mean}{\operatorname{mean}}
\newcommand{\cov}{\operatorname{cov}}
\newcommand{\dis}{\displaystyle}
% ----------------------------------------------------------------
\title{Control of surface roughness and RMS slope with patterned W}
\author{Xinyu Zhang}
% -----------------------------------------------------------------
\begin{document}
\maketitle
\tableofcontents
\section{Dynamic Model of Roughness and RMS Slope}
\subsection{2D EW equation for surface height profile}
\subsubsection{Solution via modal decomposition}
In this work, we consider a 2D Edwards-Wilkinson (EW) type equation taking the following form
\begin{equation}
\frac{\partial h}{\partial t}=w+c_2\left(\frac{\partial^2 h}{\partial x^2}+\frac{\partial^2 h}{\partial y^2}\right)+\xi(x,y,t) \label{eq:EW}
\end{equation}
where $x\in[0,L]$, $y\in[0,L]$ are the spatial coordinates, $t$ is the time, $h(x,y,t)$ is the surface height, and $\xi(x,y,t)$ is a Gaussian white noise with a zero mean and the following covariance:
\begin{equation}
\left<\xi(x,y,t)\xi(x',y',t')\right>=\sigma^2\delta(x-x')\delta(y-y')\delta(t-t').
\end{equation}
where $\delta(\cdot)$ denotes the Dirac delta function. $w$ is a patterned deposition profile with the form
\begin{equation}
\boxed{w(x,y) = w_0+A\sin\left(\frac{2k\pi}{L}x\right)},
\end{equation}
where $w_0$ is the spacial mean value of deposition rate, $A$ is the magnitude of sin wave, $k$ is an integer that responding to how many complete sin wave are there within the domain. $c_2$, $\sigma^2$ are model parameters that depends on $w_0$. The stochastic PDE of eq.~\eqref{eq:EW} is subject to the following periodic boundary conditions of the form
\begin{alignat}{2}
h(0,y,t) &= h(L,y,t)  &\qquad h(x,0,t) &= h(x,L,t)\label{eq:EW_BC1}\\
\frac{\partial h}{\partial x}(0,y,t) &= \frac{\partial h}{\partial x}(L,y,t) &\qquad \frac{\partial h}{\partial y}(x,0,t) &= \frac{\partial h}{\partial y}(x,L,t)\label{eq:EW_BC4}
\end{alignat}
and the initial condition
\begin{equation}
h(x,y,0)=h_0(x,y)
\end{equation}

To analyze the dynamics and obtain a finite-dimensional approximation of the EW equation, we first consider the eigenvalue problem of the linear operator of eq~\eqref{eq:EW} subject to the periodic boundary conditions of eq~\eqref{eq:EW_BC1}-~\eqref{eq:EW_BC4}:
\begin{gather}
\mathcal{A}\phi_{n_x,n_y}(x,y) = c_2(\frac{\partial^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2})\phi_{n_x,n_y}(x,y)=\lambda_{n_x,n_y}\phi_{n_x,n_y}(x,y),\label{eq:EW_operator}\\
\nabla^j\phi_{n_x,n_y}(0,y)=\nabla^j\phi_{n_x,n_y}(L,y),\quad j = 0,1\\
\nabla^j\phi_{n_x,n_y}(x,0)=\nabla^j\phi_{n_x,n_y}(x,L),\quad j = 0,1
\end{gather}
where $\lambda_{n_x,n_y}$ denotes an eigenvalue, $\phi_{n_x,n_y}$ denotes an eigenfunction, and $\nabla^j$, $j=0$, $1$, denotes the gradient of a given function. The solution of the eigenvalue problem is as follows:
\begin{align}
\lambda_{n_x,n_y}&= -\frac{4c_2\pi^2}{L^2}(n_x^2+n_y^2)\\
\phi_{1,n_x,n_y} &= \frac{2}{L}\sin(\frac{2n_x\pi}{L}x)\sin(\frac{2n_y\pi}{L}y)\\
\phi_{2,n_x,n_y} &= \left\{\begin{array}{ll}
                \frac{1}{L}                                                       & \textrm{$n_x=0$ and $n_y=0$}\\
                \frac{2}{L}\cos(\frac{2n_x\pi}{L}x)\cos(\frac{2n_y\pi}{L}y)       & \textrm{$n_x\neq0$ and $n_y\neq0$}\\
                \frac{\sqrt 2}{L}\cos(\frac{2n_x\pi}{L}x)\cos(\frac{2n_y\pi}{L}y) & \textrm{$n_x=0,n_y\neq0$ or $n_x\neq0,n_y=0$}
                \end{array}\right.\\
\phi_{3,n_x,n_y} &= \left\{\begin{array}{ll}
                0                                                           & \textrm{$n_x=0$}\\
                \frac{2}{L}\sin(\frac{2n_x\pi}{L}x)\cos(\frac{2n_y\pi}{L}y) & \textrm{$n_x\neq0$, $n_y\neq0$}\\
                \frac{\sqrt 2}{L}\sin(\frac{2n_x\pi}{L}x)                   & \textrm{$n_x\neq0$, $n_y=0$}
                \end{array}\right.\\
\phi_{4,n_x,n_y} &= \left\{\begin{array}{ll}
                0                                                           & \textrm{$n_y=0$}\\
                \frac{2}{L}\cos(\frac{2n_x\pi}{L}x)\sin(\frac{2n_y\pi}{L}y) & \textrm{$n_y\neq0$, $n_x\neq 0$}\\
                \frac{\sqrt 2}{L}\sin(\frac{2n_y\pi}{L}y)                   & \textrm{$n_y\neq0$, $n_x=0$}
                \end{array}\right.
\end{align}
The solution of the EW equation of eq~\eqref{eq:EW} can be expanded in an infinite series in terms of the eigenfunctions of the operator of eq~\eqref{eq:EW_operator} as follows:
\begin{equation}
h(x,y,t) = \sum_{n_x=0}^{+\infty}\sum_{n_y=0}^{+\infty}\phi_{1,n_x,n_y}z_{1,n_x,n_y}+\phi_{2,n_x,n_y}z_{2,n_x,n_y}+\phi_{3,n_x,n_y}z_{3,n_x,n_y}+\phi_{4,n_x,n_y}z_{4,n_x,n_y},\label{eq:EW_solution}
\end{equation}
where $z_{1,n_x,n_y}$, $z_{2,n_x,n_y}$, $z_{3,n_x,n_y}$, and $z_{4,n_x,n_y}$ are time-varying coefficients.

Substituting the above expansion for the solution, $h(x,t)$, into eq~\eqref{eq:EW} and taking the inner product with the adjoint eigenfunctions, the following system of infinite stochastic linear ordinary differential equations (ODEs) for the temporal evolution of the time-varying coefficients is obtained:
\begin{gather}
\frac{d z_{2,0,0}}{dt} = w_{2,0,0}+\xi_{2,0,0}(t),\label{eq:mode0_ODE}\\
\frac{d z_{p,n_x,n_y}}{dt} = w_{p,n_x,n_y}+\lambda_{n_x,n_y}z_{p,n_x,n_y}+\xi_{p,n_x,n_y}(t)\\
p = 1,2,3,4, \quad n_x,n_y = 0,1,\cdots,\infty,\quad n_x^2+n_y^2\neq0, \label{eq:mode_ODE}\notag
\end{gather}
where $\xi_{p,n_x,n_y}(t)=\iint_0^{L}\xi(x,y,t)\phi_{p,n_x,n_y}(x,y)dxdy$ is the projection of the noise $\xi(x,y,t)$ on the ODE for $z_{p,n_x,n_y}$. The noise term, $\xi_{p,n_x,n_y}$, has zero mean and covariance
\begin{equation}
\left<\xi_{p,n_x,n_y}(t)\xi_{p,n_x,n_y}(t')\right>=\sigma^2\delta(t-t').
\end{equation}
Similarly, $w_{p,n_x,n_y}$ is the projection of $w$ on the ODE for $z_{p,n_x,n_y}$, $w_{p,n_x,n_y}=\iint_{0}^{L}\phi_{p,n_x,n_y}(x,y)w(x,y)dxdy$
\begin{itemize}
% w_{1,n_x,n_y} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{If $p=1$
\begin{align}
w_{1,n_x,n_y} &= \boxed{0}
\end{align}
}
% w_{2,n_x,n_y} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{If $p=2$
    \begin{itemize}
    \item{If $n_y=0$
        \begin{itemize}
        \item{If $n_x=0$
            \begin{align}
            w_{2,0,0}     &= \iint_{0}^{L}\frac{1}{L}\left[w_0+A\sin\left(\frac{2k\pi}{L}x\right)\right]dxdy\notag\\
                          &= \boxed{w_0L+\frac{AL}{2k\pi}\left[1-\cos(2k\pi)\right]}\\
                          &= \boxed{w_0L}
            \end{align}
        }
        \item{If $n_x\neq0$
            \begin{align}
            w_{2,n_x,0} &= \iint_{0}^{L}\frac{\sqrt{2}}{L}\cos(\frac{2n_x\pi}{L}x)\left[w_0+A\sin\left(\frac{2k\pi}{L}x\right)\right]dxdy\notag\\
                        &= \sqrt{2}A\int_{0}^{L}\cos(\frac{2n_x\pi}{L}x)\sin(\frac{2k\pi}{L}x)dx\notag\\
                        &= \sqrt{2}A\int_{0}^{L}\frac{1}{2}\left(\sin(\frac{2(n_x+k)\pi}{L}x)-\sin(\frac{2(n_x-k)\pi}{L}x)\right)dx\notag
            \end{align}
            \begin{itemize}
            \item{If $n_x=k$
                \begin{align}
                w_{2,n_x,0}   &= \sqrt{2}A\int_{0}^{L}\frac{1}{2}\left(\sin(\frac{4n_x\pi}{L}x)\right)dx = \boxed{0}
                \end{align}
            }
            \item{If $n_x\neq k$ and $n_x\neq0$
                \begin{align}
                w_{2,n_x,0} &= \frac{\sqrt{2}AL}{4\pi}\left(-\frac{1}{n_x+k}\left(\cos(2(n_x+k)\pi)-1\right)+\frac{1}{n_x-k}\left(\cos(2(n_x-k)\pi)-1\right)\right)\notag\\
                            &= \boxed{\frac{\sqrt{2}ALk}{2\pi(n_x^2-k^2)}\left[\cos(2k\pi)-1\right]} = \boxed{0}
                \end{align}
            }
            \end{itemize}
        }
        \end{itemize}
    }
    \item{If $n_y\neq 0$
        \begin{align}
        w_{2,n_x,n_y} = \boxed{0}
        \end{align}
    }
    \end{itemize}
}
% w_{3,n_x,n_y} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{If $p=3$
\begin{itemize}
\item{If $n_y\neq 0$
    \begin{itemize}
    \item{If $n_x=0$
        \begin{align}
        w_{3,0,0}     &= \boxed{0}
        \end{align}
    }
    \item{If $n_x\neq 0$
        \begin{align}
        w_{3,n_x,0}   &= \iint_{0}^{L}\frac{\sqrt{2}}{L}\sin(\frac{2n_x\pi}{L}x)\left[w_0+A\sin(\frac{2k\pi}{L}x)\right]dxdy\notag\\
                      &= \sqrt{2}A\int_{0}^{L}\sin\left(\frac{2k\pi}{L}x\right)\sin\left(\frac{2m\pi}{L}x\right)dx\\
        \intertext{If $n_x=k$}
        w_{3,n_x,0}   &= \sqrt{2}A\int_{0}^{L}\sin^2\left(\frac{2n_x\pi}{L}x\right)dx = \sqrt{2}w_0\int_{0}^{L}\frac{1-\cos(4n_x\pi x/L)}{2}dx\notag = \boxed{\frac{\sqrt{2}AL}{2}}\\
        \intertext{If $n_x\neq k$}
        w_{3,n_x,0}   &= \frac{\sqrt{2}A}{2}\int_{0}^{L}\left[-\cos\left(\frac{2(k+n_x)\pi x}{L}\right)+\cos\left(\frac{2(k-n_x)\pi x}{L}\right)\right]dx\notag\\
                      &= \frac{\sqrt{2}AL}{4\pi}\left[-\frac{1}{k+n_x}\sin(2(k+n_x)\pi)+\frac{1}{k-n_x}\sin(2(k-n_x)\pi)\right]\notag\\
                      &= \boxed{\frac{\sqrt{2}ALn_x}{2\pi(k^2-n_x^2)}\sin(2k\pi)} = \boxed{0}
        \end{align}
    }
    \end{itemize}
}
\item{If $n_y=0$
    \begin{align}
    w_{3,n_x,n_y} &= \boxed{0}
    \end{align}
}
\end{itemize}
}
% w_{4,n_x,n_y} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{If $p=4$,
\begin{align}
w_{4,m,n}     &= \boxed{0}
\end{align}
}
\end{itemize}
In summary, most of the projection $w_{p,m,n}$ is 0 with two exceptions
\begin{align}
w_{1,n_x,n_y} &= 0\\
w_{2,n_x,n_y} &= \begin{cases}
                 w_0L & \textrm{$n_x=0,n_y=0$}\\
                 0    & \textrm{otherwise}
                 \end{cases}\\
w_{3,n_x,n_y} &= \begin{cases}
                 \frac{\sqrt{2}AL}{2} & \textrm{$n_x=k,n_y=0$}\\
                 0                    & \textrm{otherwise}
                 \end{cases}\\
w_{4,n_x,n_y} &= 0
\end{align}

\subsubsection{Numerical solution}
Because the ODEs eq.~\eqref{eq:mode0_ODE}, \eqref{eq:mode_ODE} are linear, we use their analytical solutions in the integration steps
\begin{align}
z_{2,0,0}(t+\Delta t)=z_{2,0,0}(t)+w_{2,0,0}\Delta t + \sigma_{2,0,0} \sqrt{\Delta t}\xi(t)\label{eq:mode0_numerical_solution}
\end{align}
\begin{align}
z_{p,n_x,n_y}(t+\Delta t) = e^{\lambda \Delta t}z_{p,n_x,n_y}(t)+\frac{e^{\lambda \Delta t-1}}{\lambda}w_p(t)+\sigma_{p,n_x,n_y} \sqrt{\frac{e^{2\lambda \Delta t}-1}{2\lambda}}\xi(t)\label{eq:mode_numerical_solution}
\end{align}
where $\xi(t)\sim N(0,1)$ is a gaussian random variable. $\xi(t)$ and $\xi(t+\Delta t)$ are independent. $\xi(t)$ is generated from gaussian random number generator.

%\begin{rem}
%One variance of equation~\eqref{eq:mode_numerical_solution} is
%\begin{align}
%z_{p,n_x,n_y}(t+\Delta t) = e^{\lambda \Delta t}z_{p,n_x,n_y}(t)+\frac{e^{\lambda \Delta t-1}}{\lambda}w_p(t)+\sigma_{p,n_x,n_y} \sqrt{\frac{e^{2\lambda \Delta t}-1}{2\lambda}}\xi(t).
%\end{align}
%This form has larger error and requires $\Delta t$ to be small enough.
%\end{rem}

\subsubsection{Dynamics of the mean and variance of states}
The temporal evolution of the variance of mode $z_{p,n_x,n_y}$ can be obtained from the solution of the linear ODEs of eqs~\eqref{eq:mode0_ODE} and \eqref{eq:mode_ODE} as follows:
\begin{align}
\left<z_{2,0,0}(t)\right> &= w_{2,0,0}(t-t_0)\label{eq:mode0_mean_dyn}\\
\var(z_{2,0,0}(t)        &= \sigma^2(t-t_0)\label{eq:mode0_cov_dyn}\\
\left<z(t)\right>        &= e^{\lambda(t-t_0)}\left<z(t_0)\right>+\frac{w_p}{\lambda}(e^{\lambda(t-t_0)}-1)\label{eq:mode_mean_dyn}\\
\var(z(t))               &= e^{2\lambda(t-t_0)}\var(z(t_0))+\sigma^2\frac{e^{2\lambda(t-t_0)}-1}{2\lambda}\label{eq:mode_cov_dyn}
\end{align}
where $z(t)=z_{p,n_x,n_y}(t)$ and $w_p=w_{p,n_x,n_y}$ for $n_x^2+n_y^2\neq0$.
\begin{rem}
The equivalence between eq.~\eqref{eq:mode0_mean_dyn} -- ~\eqref{eq:mode_cov_dyn} and eq.~\eqref{eq:mode0_numerical_solution} -- ~\eqref{eq:mode_numerical_solution} can be verified.
%\begin{equation}
%\begin{split}
%\left<z(t+\Delta t)^2\right> = \left<\left(e^{\lambda \Delta t}z(t)+\frac{e^{\lambda \Delta t-1}}{\lambda}w(t)+\sigma \sqrt{\frac{e^{2\lambda \Delta t}-1}{2\lambda}}\xi(t)\right)^2\right>\\
%&=\left<e^{2\lambda \Delta t}z^2(t)+\frac{e^{2\lambda \Delta t-2e^{\lambda \Delta}+1}}{\lambda^2}w^2(t)+\sigma^2\frac{e^{2\lambda \Delta t}-1}{2\lambda}\xi^2(t)\right.\right.\\
%&+2e^{\lambda \Delta t}z(t)\frac{e^{\lambda \Delta t-1}}{\lambda}w(t)\\
%&+2e^{\lambda \Delta t}z(t)\sigma \sqrt{\frac{e^{2\lambda \Delta t}-1}{2\lambda}}\xi(t)\\
%&\left.\left.+2\frac{e^{\lambda \Delta t-1}}{\lambda}w(t)\sigma \sqrt{\frac{e^{2\lambda \Delta t}-1}{2\lambda}}\xi(t)\right>
%\end{split}
%\end{equation}
\end{rem}
%==========================================================================================================
\subsection{Surface Roughness}
\subsubsection{Basic definition of roughness and calculation of roughness in terms of modes}
Surface roughness of the thin film is defined as the standard deviation of the surface height profile from its average height
\begin{equation}
r(t) = \sqrt{\frac{1}{L^2}\iint_{0}^{L}\left[h(x,y,t)-\bar{h}(x,y,t)\right]^2dxdy}\label{eq:r_def}
\end{equation}
where $\bar{h}(t) = \frac{1}{L^2}\iint_{0}^{L}h(x,y,t)dxdy$ is the average surface height. According to Eq.\eqref{eq:EW_solution}, we have
\begin{equation}
\bar{h}(t) = \frac{1}{L^2}\iint_{0}^{L}\phi_{2,0,0}z_{2,0,0}dxdy=\frac{1}{L}z_{2,0,0}.
\end{equation}
$\left<r^2(t)\right>$ can be rewritten as
\begin{align}
\left<r^2(t)\right>&=\left<\frac{1}{L^2}\iint_{0}^{L}\left[h(x,y,t)-\bar{h}(x,y,t)\right]dxdy\right>\label{eq:r2_mean_def}
\end{align}
where
\begin{equation}
h(x,y,t)-\bar{h}(x,y,t)=\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq0}}^{\infty}\sum_{p=1}^{4}{\phi_{p,n_x,n_y}z_{p,n_x,n_y}}
\end{equation}
\begin{align}
\begin{split}
\left<r^2(t)\right>&=\left<\frac{1}{L^2}\iint_{0}^{L}\left[\sum_{p=1}^{4}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq0}}^{\infty}z_{p,n_x,n_y}\phi_{p,n_x,n_y}(x,y)\right]^2dxdy\right>\\
                   &=\frac{1}{L^2}\iint_{0}^{L}\sum_{p_1=1}^{4}\sum_{p_2=1}^{4}\sum_{\substack{n_{1x},n_{1y}=0\\n_{1x}^2+n_{1y}^2\neq 0}}^{\infty}\sum_{\substack{n_{2x},n_{2y}=0\\n_{2x}^2+n_{2y}^2\neq 0}}^{\infty}\left<z_{p_1,n_{1x},n_{1y}}z_{p_2,n_{2x},n_{2y}}\right>\phi_{p_1,n_{1x},n_{1y}}\phi_{p_2,n_{2x},n_{2y}}dxdy\\
                   &=\frac{1}{L^2}\sum_{p_1=1}^{4}\sum_{p_2=1}^{4}\sum_{\substack{n_{1x},n_{1y}=0\\n_{1x}^2+n_{1y}^2\neq 0}}^{\infty}\sum_{\substack{n_{2x},n_{2y}=0\\n_{2x}^2+n_{2y}^2\neq 0}}^{\infty}\left<z_{p_1,n_{1x},n_{1y}}z_{p_2,n_{2x},n_{2y}}\right>\iint_{0}^{L}\phi_{p_1,n_{1x},n_{1y}}\phi_{p_2,n_{2x},n_{2y}}dxdy\\
                   &=\frac{1}{L^2}\sum_{p=1}^{4}\sum_{\substack{n_{x},n_{y}=0\\n_{x}^2+n_{y}^2\neq 0}}^{\infty}\left<z^2_{p,n_{x},n_{y}}\right>\iint_{0}^{L}\phi^2_{p,n_{x},n_{y}}dxdy\\
                   &=\boxed{\frac{1}{L^2}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq0}}^{\infty}\sum_{p=1}^{4}\left<z^2_{p,n_x,n_y}\right>}
\end{split}\label{eq:r2_mean}
\end{align}
where
\begin{equation}
\boxed{
\left<z^2_{p,n_x,n_y}\right>=\var(z_{p,n_x,n_y})+\left<z_{p,n_x,n_y}\right>^2}.
\end{equation}

%==========================================================================================================
\subsubsection{A straight-forward method to calculate aggregated roughness in 1D case}
In 1D case, the aggregated surface is defined by
\begin{equation}
\tilde{h}(\tilde{i}) = \frac{1}{l_a}\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}h(i),\quad i=0,1,\dots,\lfloor l/l_a \rfloor \label{eq:h_agg_def_1D},
\end{equation}
where $l_a$ is aggregation length (site), $l$ is the total length of the domain (site). To calculate roughness according to its definition, we can expend eq.~\eqref{eq:h_agg_def_1D} as the sum of infinite series
\begin{align}
\tilde{h}(\tilde{i}) &= \frac{1}{l_a}\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\sum_{n_x=0}^{+\infty}\left(\phi_{n_x}(i)\alpha_{n_x}(t)+\psi_{n_x}(i)\beta_{n_x}(t)\right)\label{eq:h_agg_def_expended_1D}
\end{align}
Roughness is defined the same way as
\begin{align}
\begin{split}
\displaystyle\tilde{r}^2 &= \frac{1}{(l/l_a)}
\sum_{\tilde{i}=0}^{l/l_a}\left(\tilde{h}(\tilde{i})-\bar{h}\right)^2
\end{split}
\end{align}

To have an estimation of the computational complexity of the above problem, consider the following example: $L = 8000$ nm, $l = 8000/0.2 = 40000$ ($0.2$ nm per point on the ``high resolution'' surface), $n_x = 0,1,\dots,50$, $i=0,1,\dots,399999$, $l_a=400$, thus $\tilde{i}= 0,\dots,99$. Then, for each $\tilde{i}$, eq~\eqref{eq:h_agg_def_expended_1D} needs $400\times50\times2=4.0\times10^4$ operations. To get the whole aggregated surface, $100\times4.0\times10^4=\boxed{4.0\times10^{6}}$ operations are needed.

\subsubsection{An improved method to calculate aggregated roughness in 1D case}
To speedup the calculation of roughness, we can use the following formula
\begin{align}
\begin{split}
\tilde{h}(\tilde{i}) &= \frac{1}{l_a}\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\sum_{n_x=0}^{+\infty}\left(\phi_{n_x}(i)\alpha_{n_x}(t)+\psi_{n_x}(i)\beta_{n_x}(t)\right)\\
             &=\sum_{n_x=0}^{+\infty}\alpha_{n_x}(t)\left(\frac{1}{l_a}\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\phi_{n_x}(i)\right)+\beta_{n_x}(t)\left(\frac{1}{l_a}\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\psi_{n_x}(i)\right)\\
             &=\sum_{n_x=0}^{+\infty}A_{n_x}(\tilde{i})\alpha_{n_x}(t)\phi_{n_x}(\tilde{i})+B_{n_x}(\tilde{i})\beta_{n_x}(t)\psi_{n_x}(\tilde{i})
\end{split}\label{eq:h_1D_accelerate}
\end{align}

\begin{rem}
With eq~\eqref{eq:h_1D_accelerate}, we cannot calculate roughness like eq~\eqref{eq:r2_mean}
\begin{equation}
\tilde{r}^2(t)\neq \sum_{n_x=0}^{+\infty}\left<A_{n_x}(\tilde{i})\alpha_{n_x}(t)\right>^2+\left<B_{n_x}(\tilde{i})\beta_{n_x}(t)\right>^2
\end{equation}
because I cannot prove the orthogonality of $A_{n_x}(\tilde{i})\phi_{n_x}(\tilde{i})$.
\end{rem}
%==========================================================================================================
\subsubsection{A straight-forward method to calculate aggregated roughness in 2D case}
The aggregated surface is defined by
\begin{equation}
\tilde{h}(\tilde{i},\tilde{j}) = \frac{1}{l_a^2}\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\sum_{j=\tilde{j}*l_a}^{(\tilde{j}+1)*l_a-1}h(i,j),\quad i,j=0,1,\dots,\lfloor l/l_a \rfloor \label{eq:h_agg_def}
\end{equation}

First we show that direct calculation aggregated roughness based on surface profile might be impractical. Expend eq.~\eqref{eq:h_agg_def} using eq.~\eqref{eq:EW_solution}
\begin{align}
\tilde{h}(\tilde{i},\tilde{j}) &= \frac{1}{l_a^2}\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\sum_{j=\tilde{j}*l_a}^{(\tilde{j}+1)*l_a-1}\left(\sum_{n_x=0}^{+\infty}\sum_{n_y=0}^{+\infty}\sum_{p=1}^{4}\phi_{p,n_x,n_y}(x_i,y_j)z_{p,n_x,n_y}(t)\right)\label{eq:h_agg_def_expended}
\end{align}
Roughness is defined the same way as eq.~\eqref{eq:r2_mean_def}
\begin{align}
\begin{split}
\displaystyle\tilde{r}^2 &= \frac{1}{(l/l_a)^2}
\sum_{\tilde{i}=0}^{l/l_a}\sum_{\tilde{j}=0}^{l/l_a}\left(\tilde{h}(\tilde{i},\tilde{j})-\bar{h}\right)^2
\end{split}
\end{align}

To have an estimation of the computational complexity of the above problem, consider the following example: $L = 8000$ nm, $l = 8000/0.2 = 40000$ ($0.2$ nm per point on the ``high resolution'' surface), $n_x,n_y = 0,1,\dots,50$, $i,j=0,1,\dots,399999$, $l_a=400$ ($100\times100$ points on the aggregated surface), thus $\tilde{i},\tilde{j}=0,\dots, 99$.

For each $(\tilde{i},\tilde{j})$, eq~\eqref{eq:h_agg_def_expended} needs $400\times400\times50\times50\times4=1.6\times10^9$ operations. To get the whole aggregated surface, $100\times100\times1.6\times10^9=\boxed{1.6\times10^{13}}$ operations are needed. This calculation load is heavy for open simulations and it might be impractical to to be used as a model for MPC controller. Assume we use finite difference to approximate the gradient, then each gradient calculation needs at least 3 function value calculation. Each iteration needs at least one gradient calculation and one function value calculation (similar to Newton's method when solving nonlinear equation). Assume each optimization step needs 10 iterations, then we need $10\times (1+3)\times 1.6\times10^{13}=6.4\times10^{14}$ operations to solve one step of MPC.
%==========================================================================================================
\subsubsection{An Improved method to calculate aggregated surface and rms slope}
One possible solution of this problem is to calculate roughness directly from mode without constructing the detailed surface profile.
\begin{align}
\tilde{h}(\tilde{i},\tilde{j}) &= \frac{1}{l_a^2}\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\sum_{j=\tilde{j}*l_a}^{(\tilde{j}+1)*l_a-1}\left(\sum_{n_x=0}^{+\infty}\sum_{n_y=0}^{+\infty}\sum_{p=1}^{4}\phi_{p,n_x,n_y}(x_i,y_j)z_{p,n_x,n_y}(t)\right)\notag\\
             &=\sum_{n_x=0}^{+\infty}\sum_{n_y=0}^{+\infty}\sum_{p=1}^{4}z_{p,n_x,n_y}(t)\left(\frac{1}{l_a^2}\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\sum_{j=\tilde{j}*l_a}^{(\tilde{j}+1)*l_a-1}\phi_{p,n_x,n_y}(x_i,y_j)\right)\notag\\
             &=\sum_{n_x=0}^{+\infty}\sum_{n_y=0}^{+\infty}\sum_{p=1}^{4}A_{p,n_x,n_y}(\tilde{i},\tilde{j})z_{p,n_x,n_y}(t)\phi_{p,n_x,n_y}(x_{\tilde{i}},y_{\tilde{j}})
\end{align}
where
\begin{align}
A_{p,n_x,n_y}(\tilde{i},\tilde{j}) &= \left(\frac{1}{l_a^2\phi_{p,n_x,n_y}(x_{\tilde{i}},y_{\tilde{j}})}\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\sum_{j=\tilde{j}*l_a}^{(\tilde{j}+1)*l_a-1}\phi_{p,n_x,n_y}(x_i,y_j)\right)
\end{align}

If $A_{p,n_x,n_y}(\tilde{i},\tilde{j})$ is independent of $(\tilde{i},\tilde{j})$, then eq.~\eqref{}
\begin{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{For $p=1$
\begin{align}
A_{1,n_x,n_y}(\tilde{i},\tilde{j}) &= \left(\frac{1}{l_a^2\phi_{p,n_x,n_y}(x_{\tilde{i}},y_{\tilde{j}})}\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\sum_{j=\tilde{j}*l_a}^{(\tilde{j}+1)*l_a-1}\frac{2}{L}\sin(\frac{2n_x\pi}{l}i)\sin(\frac{2n_y\pi}{l}j)\right)\notag\\
                                   &=\frac{1}{l_a^2\sin(\frac{2n_x\pi x_{\tilde{i}}}{L})\sin(\frac{2n_y\pi y_{\tilde{j}}}{L})}\left(\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\sin(\frac{2n_x\pi}{l}i)\right)\left(\sum_{j=\tilde{j}*l_a}^{(\tilde{j}+1)*l_a-1}\sin(\frac{2n_y\pi}{l}j)\right)\notag\\
                                   &=\left(\frac{1}{l_a\sin(\frac{2n_x\pi x_{\tilde{i}}}{L})}\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\sin(\frac{2n_x\pi}{l}i)\right)
                                   \left(\frac{1}{l_a\sin(\frac{2n_y\pi y_{\tilde{j}}}{L})}\sum_{j=\tilde{j}*l_a}^{(\tilde{j}+1)*l_a-1}\sin(\frac{2n_y\pi}{l}j)\right)\label{eq:A1mn}
\end{align}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{For $p=2$
\begin{align}
&A_{1,n_x,n_y}(\tilde{i},\tilde{j}) =\notag\\
&\begin{cases}
1 & n_x=0 \text{ and } n_y=0\\
\displaystyle\frac{1}{l_a^2\cos(\frac{2n_x\pi x_{\tilde{i}}}{L})\cos(\frac{2n_y\pi y_{\tilde{j}}}{L})}\left(\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\cos(\frac{2n_x\pi}{l}i)\right)\left(\sum_{j=\tilde{j}*l_a}^{(\tilde{j}+1)*l_a-1}\cos(\frac{2n_y\pi}{l}j)\right) & n_x\neq 0 \text{ and } n_y\neq 0\\
\displaystyle\frac{1}{l_a^2\cos(\frac{2n_x\pi x_{\tilde{i}}}{L})\cos(\frac{2n_y\pi y_{\tilde{j}}}{L})}\left(\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\cos(\frac{2n_x\pi}{l}i)\right)\left(\sum_{j=\tilde{j}*l_a}^{(\tilde{j}+1)*l_a-1}\cos(\frac{2n_y\pi}{l}j)\right) & n_x = 0,n_y\neq 0 \text{ or } n_x \neq 0,n_y = 0
\end{cases}
\end{align}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{For $p=3$
\begin{align}
&A_{1,n_x,n_y}(\tilde{i},\tilde{j}) =\notag\\
&\begin{cases}
0 & n_x=0 \text{ and } n_x=0\\
\displaystyle\frac{1}{l_a^2\sin(\frac{2n_x\pi x_{\tilde{i}}}{L})\cos(\frac{2n_y\pi y_{\tilde{j}}}{L})}\left(\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\sin(\frac{2n_x\pi}{l}i)\right)\left(\sum_{j=\tilde{j}*l_a}^{(\tilde{j}+1)*l_a-1}\cos(\frac{2n_y\pi}{l}j)\right) & n_x\neq 0, n_y\neq 0\\
\displaystyle\frac{1}{l_a\sin(\frac{2n_x\pi x_{\tilde{i}}}{L})}\left(\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\sin(\frac{2n_x\pi}{l}i)\right) & n_x \neq 0,n_y = 0
\end{cases}
\end{align}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{For $p=4$
\begin{align}
&A_{1,n_x,n_y}(\tilde{i},\tilde{j}) =\notag\\
&\begin{cases}
0 & n_x=0 \text{ and } n_y=0\\
\displaystyle\frac{1}{l_a^2\cos(\frac{2n_x\pi x_{\tilde{i}}}{L})\sin(\frac{2n_y\pi y_{\tilde{j}}}{L})}\left(\sum_{i=\tilde{i}*l_a}^{(\tilde{i}+1)*l_a-1}\cos(\frac{2n_x\pi}{l}i)\right)\left(\sum_{j=\tilde{j}*l_a}^{(\tilde{j}+1)*l_a-1}\sin(\frac{2n_y\pi}{l}j)\right) & n_x\neq 0, n_y\neq 0\\
\displaystyle\frac{\sqrt{1}}{l_a^2\sin(\frac{2n_y\pi y_{\tilde{j}}}{L})}\left(\sum_{j=\tilde{j}*l_a}^{(\tilde{j}+1)*l_a-1}\frac{2}{L}\sin(\frac{2n_y\pi}{l}j)\right) & n_x = 0,n_y\neq 0
\end{cases}
\end{align}
}
\end{itemize}

Then, the eigenvalue spectrum of $\tilde{h}$ is $\tilde{z}_{p,n_x,n_y}=A_{p,n_x,n_y}z_{p,n_x,n_y}$. Roughness and slope can be calculated according to eq.~\eqref{eq:r2_mean} and ~\eqref{eq:m2_mean}

Eq.~\eqref{eq:A1mn} can be further simplified
\begin{align}
&\displaystyle\frac{1}{\sin(\frac{2 n_x \pi \tilde{i}l_a}{l})}\sum_{i=\tilde{i}l_a}^{\tilde{i}l_a+l_a-1}\sin(\frac{2n_x\pi}{l}i)=\frac{1}{\sin(\frac{2 n_x \pi \tilde{i}l_a}{l})}\sum_{I=0}^{l_a-1}\sin(\frac{2n_x\pi}{l}(\tilde{i}l_a+I))\notag\\
=&\frac{1}{\sin(\frac{2 n_x \pi \tilde{i}l_a}{l})}\left(\sum_{I=0}^{l_a-1}
\sin(\frac{2n_x\pi}{l}\tilde{i}l_a)\cos(\frac{2n_x\pi}{l}I)+
\cos(\frac{2n_x\pi}{l}\tilde{i}l_a)\sin(\frac{2n_x\pi}{l}I)\right)\notag\\
=&\sum_{I=0}^{l_a-1}
\cos(\frac{2n_x\pi}{l}I)+
\left(\sum_{I=0}^{l_a-1}\sin(\frac{2n_x\pi}{l}I)\right)/\tan(\frac{2n_x\pi}{l}\tilde{i}l_a)\notag\\
&\tilde{i} = 0,1,\dots,\lfloor l/la \rfloor \notag
\end{align}

The roughness is defined similarly as the original roughness,
\begin{equation}
\begin{split}
\tilde{r}(t)&= \sqrt{\frac{1}{L^2}\sum_{\tilde{i}=0}^{L}\sum_{\tilde{j}=0}^{L}\left[\tilde{h}(x,y,t)-\bar{h}(x,y,t)\right]^2}\\
    &= \sqrt{\frac{1}{L^2}\sum_{\tilde{i}=0}^{L}\sum_{\tilde{j}=0}^{L}\left[\tilde{h}(x,y,t)-\bar{h}(x,y,t)\right]^2}
\end{split}
\end{equation}

%\subsection{Roughness and slope of re-sampled surface}
%\section{Detailed implementation}
%\section{Result validation}
%==========================================================================================================
\subsection{RMS slope}
The rms slope is defined as the root--mean--square of the slope of the surface height
\begin{gather}
\begin{split}
m(t) &= \sqrt{\frac{1}{L^2}\iint_0^{L}\left(\frac{\partial h}{\partial x}\right)^2dxdy}\\
     &= \sqrt{\frac{1}{L^2}\sum_{i=0}^{l-1}\sum_{j=0}^{l-1}\left(\frac{\hat{h}(i+1,j)-\hat{h}(i,j)}{\Delta x}\right)^2\frac{L^2}{l^2}}
\end{split}
\end{gather}
\begin{gather}
\begin{split}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\left<m^2(t)\right>&=\left<\frac{1}{L^2}\sum_{i=0}^{l-1}\sum_{j=0}^{l-1}\left(\frac{\hat{h}(i+1,j)-\hat{h}(i,j)}{\Delta x}\right)^2\frac{L^2}{l^2}\right>\\
                   &=\left<\frac{1}{L^2}\sum_{i=0}^{l-1}\sum_{j=0}^{l-1}\left(\sum_{p=1}^{4}\sum_{n_x,n_y=0}z_{p,n_x,n_y}\left(\phi_{p,n_x,n_y}(i+1,j)-\phi_{p,n_x,n_y}(i,j)\right)\right)^2\right>\\
                   &=\left<\frac{1}{L^2}\sum_{i=0}^{l-1}\sum_{j=0}^{l-1}\sum_{p_1=1}^{4}\sum_{\substack{n_{1x}=0\\n_{1y}=0}}\sum_{p_2=1}^{4}\sum_{\substack{n_{2x}=0\\n_{2y}=0}}z_{p_1,n_{1x},n_{1y}}z_{p_2,n_{2x},n_{2y}}\Delta \phi_{p_1,n_{1x},n_{1y}}(i,j)\Delta\phi_{p_2,n_{2x},n_{2y}}(i,j)\right>\\
                   &=\frac{1}{L^2}\sum_{p_1=1}^{4}\sum_{\substack{n_{1x}=0\\n_{1y}=0}}\sum_{p_2=1}^{4}\sum_{\substack{n_{2x}=0\\n_{2y}=0}}\left<z_{p_1,n_{1x},n_{1y}}z_{p_2,n_{2x},n_{2y}}\right>\left(\sum_{i=0}^{l-1}\sum_{j=0}^{l-1}\Delta \phi_{p_1,n_{1x},n_{1y}}(i,j)\Delta\phi_{p_2,n_{2x},n_{2y}}(i,j)\right)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\left<m^2(t)\right>&=\left<\frac{1}{L^2}\iint_{0}^{L}\left[\sum_{p=1}^{4}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq0}}^{\infty}z_{p,n_x,n_y}\frac{\partial \phi_{p,n_x,n_y}(x,y)}{\partial x}\right]^2dxdy\right>\\
%                   &=\frac{1}{L^2}\iint_{0}^{L}\sum_{p_1=1}^{4}\sum_{p_2=1}^{4}\sum_{\substack{n_{1x},n_{1y}=0\\n_{1x}^2+n_{1y}^2\neq 0}}^{\infty}\sum_{\substack{n_{2x},n_{2y}=0\\n_{2x}^2+n_{2y}^2\neq 0}}^{\infty}\left<z_{p_1,n_{1x},n_{1y}}z_{p_2,n_{2x},n_{2y}}\right>\frac{\partial \phi_{p_1,n_{1x},n_{1y}}}{\partial x}\frac{\partial \phi_{p_2,n_{2x},n_{2y}}}{\partial x}dxdy\\
%                   &=\frac{1}{L^2}\sum_{p_1=1}^{4}\sum_{p_2=1}^{4}\sum_{\substack{n_{1x},n_{1y}=0\\n_{1x}^2+n_{1y}^2\neq 0}}^{\infty}\sum_{\substack{n_{2x},n_{2y}=0\\n_{2x}^2+n_{2y}^2\neq 0}}^{\infty}\left<z_{p_1,n_{1x},n_{1y}}z_{p_2,n_{2x},n_{2y}}\right>\iint_{0}^{L}\frac{\partial \phi_{p_1,n_{1x},n_{1y}}}{\partial x}\frac{\partial \phi_{p_2,n_{2x},n_{2y}}}{\partial x}dxdy\\
%                   &=\frac{1}{L^2}\sum_{p=1}^{4}\sum_{\substack{n_{x},n_{y}=0\\n_{x}^2+n_{y}^2\neq 0}}^{\infty}\left<z^2_{p,n_{x},n_{y}}\right>\iint_{0}^{L}\left(\frac{\partial \phi_{p,n_{x},n_{y}}}{\partial x}\right)^2dxdy\\
%                   &=\frac{1}{L^2}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq0}}^{\infty}\sum_{p=1}^{4}\left<z^2_{p,n_x,n_y}\right>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{split}
\end{gather}
When $p_1\neq p_2$ or $n_{1x}\neq n_{2x}$ or $n_{1y}\neq n_{2y}$, $\left<z_{p_1,n_{1x},n_{1y}}z_{p_2,n_{2x},n_{2y}}\right> = 0$ because $z_{p_1,n_{1x},n_{1y}}$ and $z_{p_2,n_{2x},n_{2y}}$ are independent. Thus
\begin{equation}
\left<m^2(t)\right> = \frac{1}{L^2}\sum_{p=1}^{4}\sum_{\substack{n_{x}=0\\n_{y}=0}}\left<z^2_{p,n_{x},n_{y}}\right>\left(\sum_{i=0}^{l-1}\sum_{j=0}^{l-1}\Delta \phi^2_{p,n_{x},n_{y}}(i,j)\right)
\end{equation}

\begin{gather}
\begin{split}
 &\sum_{i=0}^{l-1}\sum_{j=0}^{l-1}\Delta \phi_{p_1,n_{1x},n_{1y}}(i,j)\Delta\phi_{p_2,n_{2x},n_{2y}}(i,j)\\
=&\sum_{i=0}^{l-1}\sum_{j=0}^{l-1}\left(\phi_{p_1,n_{1x},n_{1y}}(i+1,j)-\phi_{p_1,n_{1x},n_{1y}}(i,j)\right)\left(\phi_{p_2,n_{2x},n_{2y}}(i+1,j)-\phi_{p_2,n_{2x},n_{2y}}(i,j)\right)\\
=&\frac{4}{L^2}\left(\sum_{i=0}^{l-1}\left(\sin(\frac{2n_{1x}\pi}{l}(i+1))-\sin(\frac{2n_{1x}\pi}{l}i)\right)\left(\sin(\frac{2n_{2x}\pi}{l}(i+1))-\sin(\frac{2n_{2x}\pi}{l}i)\right)\right)\left(\sum_{j=0}^{l-1}\sin(\frac{2n_{1y}\pi}{l}j)\sin(\frac{2n_{2y}\pi}{l}j)\right)
\end{split}\\
\left<m^2(t)\right> = \sum_{p=1}^{4}\sum_{\substack{n_{x}=0\\n_{y}=0}}K_{p,n_x,n_y}\left<z_{p,n_x,n_y}^2\right>=\boxed{\sum_{p=1}^{4}\sum_{\substack{n_{x}=0\\n_{y}=0}}\frac{4}{L^2}\sin^2\left(\frac{\pi n_x}{l}\right)\left<z_{p,n_x,n_y}^2\right>}
\end{gather}
In terms of eigenfunctions
\begin{equation}
\left<m^2(t)\right>=\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{l/2}\left(K_{1,n_x,n_y}\left<z^2_{1,n_x,n_y}\right>+K_{2,n_x,n_y}\left<z^2_{2,n_x,n_y}\right>+K_{3,n_x,n_y}\left<z^2_{3,n_x,n_y}\right>+K_{4,n_x,n_y}\left<z^2_{4,n_x,n_y}\right>\right)\label{eq:m2_mean}
\end{equation}
where
\begin{equation}
K_{p,n_x,n_y} = \boxed{\frac{4}{L^2}\sin^2\left(\frac{\pi n_x}{l}\right)}
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Parameter estimation}
\subsection{Parameters for domain size 50 nm}
\subsection{Parameters for domain size 8000 nm}
\begin{equation}
c_2(l=40000) = p_1w^3+p_2w^2+p_3w+p_4
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Model Predictive Controller Design 1}
\subsection{Mathematical formulation}
MPC formulation similar to Ni Dong's and the one used in last IECR paper is used
\begin{gather}\label{eq:MPC}
\dis{\min_{T,W,A}} f(T,W,A) = q_{r^2}\left(r^2_{set}-\left<r^2_f\right>\right)^2+q_{m^2}\left(m^2_{set}-\left<m^2_f\right>\right)^2
\end{gather}
%\operatornamewithlimits{exp}_{x \in \mathbb{R}^n} J(t) = q_{r^2}\left(r^2(t_f)-r^2_{set}\right)^2+q_{m^2}\left(m^2(t_f)-m^2_{set}\right)^2\\
where
\begin{align}
\left<r^2_f\right> = \frac{1}{L^2}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{\infty}\sum_{p=1}^{4}\left<z^2_{p,n_x,n_y}(t_f)\right>,\quad
\left<m^2_f\right> = \sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}\sum_{p=1}^{4}\left(K_{p,n_x,n_y}\left<z^2_{p,n_x,n_y}(t_f)\right>\right)
\end{align}
\begin{align}
\left<z^2_{p,n_x,n_y}(t_f)\right> &=\var(z_{p,n_x,n_y}(t_f))+\left<z_{p,n_x,n_y}(t_f)\right>^2\\
\left<z_{p,n_x,n_y}(t_f)\right>   &= e^{\lambda_{n_x,n_y}(t_f-t)}\left<z_{p,n_x,n_y}(t)\right>+\frac{w_p}{\lambda_{n_x,n_y}}(e^{\lambda_{n_x,n_y}(t_f-t)}-1)\\
\var(z_{p,n_x,n_y}(t_f))     &= e^{2\lambda_{n_x,n_y}(t_f-t)}\var(z_{p,n_x,n_y}(t))+\sigma^2\frac{e^{2\lambda_{n_x,n_y}(t_f-t)}-1}{2\lambda_{n_x,n_y}}\\
\lambda_{n_x,n_y}            &= -\frac{4c_2\pi^2}{L^2}(n_x^2+n_y^2), \quad n_x^2+n_y^2\neq0
\end{align}
\begin{align}
%c(T,W)        &= W\left(1-\frac{k_w}{W^{a_w}}e^{k_BT/E_w}\right)\\
c_2(T,W)      &= \frac{k_c}{l^2W^{a_c}}e^{k_BT/E_c}\\
\sigma^2(T,W) &= \frac{\pi^2}{l^2}\left[1+e^{(a_t+k_tW)T-a_v-k_vW}\right]
\end{align}
Subject to:
\begin{gather}
T_{min}\leq T \leq T_{max}, \quad |T(t)-T(t-dt)|\leq \Delta T_{max}, \\
W_{min}\leq W \leq W_{max}, \quad |W(t)-W(t-dt)|\leq \Delta W_{max},\label{eq:MPC_LastConstraint}
\end{gather}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Implementation using the IPOPT optimizer}
In order to implement in IPOPT, the optimization problem in eq.~\eqref{eq:MPC}--\eqref{eq:MPC_LastConstraint} is rewritten as
\begin{itemize}
\item{
Cost function
\begin{gather}\label{eq:MPC}
f(T,W,A) = q_{r^2}\left(r^2_{set}-\left<r^2_f\right>\right)^2+q_{m^2}\left(m^2_{set}-\left<m^2_f\right>\right)^2
\end{gather}
%\operatornamewithlimits{exp}_{x \in \mathbb{R}^n} J(t) = q_{r^2}\left(r^2(t_f)-r^2_{set}\right)^2+q_{m^2}\left(m^2(t_f)-m^2_{set}\right)^2\\
where
\begin{align}
\left<r^2_f\right> = \frac{1}{L^2}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{\infty}\sum_{p=1}^{4}\left<z^2_{p,n_x,n_y}(t_f)\right>,\quad
\left<m^2_f\right> = \sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}\sum_{p=1}^{4}\left(K_{p,n_x,n_y}\left<z^2_{p,n_x,n_y}(t_f)\right>\right)
\end{align}
\begin{align}
\left<z^2_{p,n_x,n_y}(t_f)\right> &=\var(z_{p,n_x,n_y}(t_f))+\left<z_{p,n_x,n_y}(t_f)\right>^2\\
\left<z_{p,n_x,n_y}(t_f)\right>   &= e^{\lambda_{n_x,n_y}(t_f-t)}\left<z_{p,n_x,n_y}(t)\right>+\frac{w_p}{\lambda_{n_x,n_y}}(e^{\lambda_{n_x,n_y}(t_f-t)}-1)\\
\var(z_{p,n_x,n_y}(t_f))     &= e^{2\lambda_{n_x,n_y}(t_f-t)}\var(z_{p,n_x,n_y}(t))+\sigma^2\frac{e^{2\lambda_{n_x,n_y}(t_f-t)}-1}{2\lambda_{n_x,n_y}}\\
\lambda_{n_x,n_y}            &= -\frac{4c_2\pi^2}{L^2}(n_x^2+n_y^2), \quad n_x^2+n_y^2\neq0
\end{align}
\begin{align}
%c(T,W)        &= W\left(1-\frac{k_w}{W^{a_w}}e^{k_BT/E_w}\right)\\
c_2(T,W)      &= \frac{k_c}{l^2W^{a_c}}e^{k_BT/E_c}\\
\sigma^2(T,W) &= \frac{\pi^2}{l^2}\left[1+e^{(a_t+k_tW)T-a_v-k_vW}\right]
\end{align}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Gradient of cost function
\begin{gather}
\nabla f = \left[\frac{\partial f}{\partial T},\frac{\partial f}{\partial W},\frac{\partial f}{\partial A}\right]
\end{gather}
%%%%%%%%%%%%%%%%%%% \frac{\partial f}{\partial T} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
$\frac{\partial f}{\partial T}$ is given by
\begin{gather}
\frac{\partial f}{\partial T} = -2q_{r^2}\left(r^2_{set}-r^2_f\right)\frac{\partial r^2_f}{\partial T}-2q_{m^2}\left(m^2_{set}-m^2_f\right)\frac{\partial m^2_f}{\partial T}
\end{gather}
\begin{gather}
\frac{\partial r^2_f}{\partial T} = \frac{1}{L^2}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{\infty}\sum_{p=1}^{4}\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial T}\\
\frac{\partial m^2_f}{\partial T} = \sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{\infty}\sum_{p=1}^{4}\left(K_{p,n_x,n_y}\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial T}\right)
\end{gather}
\begin{gather}
\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial T}=\frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial T}+2\left<z_{p,n_x,n_y}\right>\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial T}
\end{gather}
\begin{gather}
\frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial T} = \frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial \lambda_{n_x,n_y}}\frac{\partial \lambda_{n_x,n_y}}{\partial T}+\frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial \sigma^2}\frac{\partial \sigma^2}{\partial T}
\end{gather}
\begin{gather}
\begin{split}
\frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial \lambda_{n_x,n_y}} =& 2e^{2\lambda_{n_x,n_y}(t_f-t)}\var(z_{p,n_x,n_y}(t))(t_f-t)\\
                                                                     &+\sigma^2\frac{4\lambda_{n_x,n_y}e^{2\lambda_{n_x,n_y}(t_f-t)}(t_f-t)-2\left(e^{2\lambda_{n_x,n_y}(t_f-t)}-1\right)}{4\lambda_{n_x,n_y}^2}
\end{split}\\
\frac{\partial \lambda_{n_x,n_y}}{\partial T}=-\frac{4\pi^2}{L^2}(n_x^2+n_y^2)\frac{\partial c_2}{\partial T}\\
\frac{\partial c_2}{\partial T} = \frac{k_c}{l^2W^{a_c}}e^{k_BT/E_c}\frac{k_B}{E_c}
\end{gather}
\begin{gather}
\frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial \sigma^2} = \frac{e^{2\lambda_{n_x,n_y}(t_f-t)}-1}{2\lambda_{n_x,n_y}}\\
\frac{\partial \sigma^2}{\partial T} = \frac{\pi^2}{l^2}e^{(a_t+k_tW)T-a_v-k_vW}(a_t+k_tW)
\end{gather}
\begin{gather}
\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial T}=\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial \lambda_{n_x,n_y}}\frac{\partial \lambda_{n_x,n_y}}{\partial T}\\
\begin{split}
\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial \lambda_{n_x,n_y}} &=e^{\lambda_{n_x,n_y}(t_f-t)}\left<z_{p,n_x,n_y}(t)\right>(t_f-t)\frac{\partial \lambda_{n_x,n_y}}{\partial T}\\
&+w_p\frac{\lambda_{n_x,n_y}e^{\lambda_{n_x,n_y}(t_f-t)}(t_f-t)-(e^{\lambda_{n_x,n_y}(t_f-t)}-1)}{\lambda^2_{n_x,n_y}}
\end{split}
\end{gather}
%%%%%%%%%%%%%%%%%%% \frac{\partial f}{\partial W} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
$\frac{\partial f}{\partial W}$ is given by
\begin{gather}
\frac{\partial f}{\partial W}= -2q_{r^2}\left(r^2_{set}-r^2_f\right)\frac{\partial r^2_f}{\partial W}-2q_{m^2}\left(m^2_{set}-m^2_f\right)\frac{\partial m^2_f}{\partial W}
\end{gather}
\begin{gather}
\frac{\partial r^2_f}{\partial W} = \frac{1}{L^2}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{\infty}\sum_{p=1}^{4}\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial W}\\
\frac{\partial m^2_f}{\partial W} = \sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{\infty}\sum_{p=1}^{4}\left(K_{p,n_x,n_y}\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial W}\right)
\end{gather}
\begin{gather}
\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial W}=\frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial W}+2\left<z_{p,n_x,n_y}\right>\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial W}
\end{gather}
\begin{gather}
\frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial W} = \frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial \lambda_{n_x,n_y}}\frac{\partial \lambda_{n_x,n_y}}{\partial W}+\frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial \sigma^2}\frac{\partial \sigma^2}{\partial W}
\end{gather}
\begin{gather}
\frac{\partial \lambda_{n_x,n_y}}{\partial W}= -\frac{4\pi^2}{L^2}(n_x^2+n_y^2)\frac{\partial c_2}{\partial W}\\
\frac{\partial c_2}{\partial W} = \frac{k_c}{l^2}e^{k_BT/E_c}\frac{-a_c}{W^{a_c+1}}
\end{gather}
\begin{gather}
\frac{\partial \sigma^2}{\partial W} = \frac{\pi^2}{l^2}e^{(a_t+k_tW)T-a_v-k_vW}(k_tT-k_v)
\end{gather}
\begin{gather}
\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial W}=\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial \lambda_{n_x,n_y}}\frac{\partial \lambda_{n_x,n_y}}{\partial W}+\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial w_{p,m,n}}\frac{\partial w_{p,m,n}}{\partial W}\\
\frac{\partial w_{p,m,n}}{\partial W}=
\begin{cases}
L    & p=2,m,n=0\\
0    & \text{otherwise}
\end{cases}
\end{gather}
%%%%%%%%%%%%%%%%%%% \frac{\partial f}{\partial A} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
$\frac{\partial f}{\partial A}$ is given by
\begin{gather}
\frac{\partial f}{\partial A}= -2q_{r^2}\left(r^2_{set}-r^2_f\right)\frac{\partial r^2_f}{\partial A}-2q_{m^2}\left(m^2_{set}-m^2_f\right)\frac{\partial m^2_f}{\partial A}
\end{gather}
\begin{gather}
\frac{\partial r^2_f}{\partial A} = \frac{1}{L^2}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{\infty}\sum_{p=1}^{4}\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial A}\\
\frac{\partial m^2_f}{\partial A} = \sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{\infty}\sum_{p=1}^{4}\left(K_{p,n_x,n_y}\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial A}\right)
\end{gather}
\begin{gather}
\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial A}=2\left<z_{p,n_x,n_y}\right>\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial A}
\end{gather}
\begin{gather}
\begin{split}
\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial A}=\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial w_p}\frac{\partial w_p}{\partial A} = \frac{1}{\lambda_{n_x,n_y}}\left(e^{\lambda_{n_x,n_y}(t_f-t)}-1\right)\frac{\partial w_{p,n_x,n_y}}{\partial A}
\end{split}\\
\frac{\partial w_{p,n_x,n_y}}{\partial A} =
\begin{cases}
0                                         & \text{if $p=1$}\\
\frac{L}{2k\pi}\left[1-\cos(2k\pi)\right] & \text{if $p=2, n_x,n_y=0$}\\
0                                         & \text{if $p=2, n_x,n_y \neq 0$}\\
0                                         & \text{if $p=2, n_x=k, n_y  = 0$}\\
\frac{\sqrt{2}Lk}{2\pi(n_x^2-k^2)}\left[\cos(2k\pi)-1\right] & \text{if $p=2, n_x\neq k, n_y = 0$}\\
0                                         & \text{if $p=3, n_x,n_y = 0$}\\
0                                         & \text{if $p=3, n_x,n_y \neq 0$}\\
\frac{\sqrt{2}L}{2}                       & \text{if $p=3, n_x=k, n_y =0$}\\
\frac{\sqrt{2}Ln_x}{2\pi(k^2-n_x^2)\sin(2k\pi)} & \text{if $p=3, n_x\neq k, n_y=0$}\\
0                                         & \text{if $p=4$}
\end{cases}
\end{gather}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Inequality constraints
\begin{gather}
0\leq g(x) = W-A \leq \infty
\end{gather}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Jacobian matrix of constraints
\begin{gather}
\nabla g = \begin{bmatrix} 0 & 1 & -1 \end{bmatrix}
\end{gather}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Variable constraints
\begin{gather}
\max\{T_{min},T_{old}-\Delta T_{max}\}\leq T \leq \min\{T_{max},T_{old}+\Delta T_{max}\}\\
\max\{W_{min},W_{old}-\Delta W_{max}\}\leq W \leq \min\{W_{max},W_{old}+\Delta W_{max}\}\\
0 \leq A \leq \infty
\end{gather}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Hessian matrix of the Lagrangian
\begin{gather}
\sigma_f\nabla^2 f+\lambda_1\nabla^2 g=
\sigma_f
\begin{bmatrix}
\frac{\partial^2 f}{\partial T^2}         &                                            & \\
\frac{\partial^2 f}{\partial W\partial T} & \frac{\partial^2 f}{\partial W^2}          & \\
\frac{\partial^2 f}{\partial A\partial T} & \frac{\partial^2 f}{\partial A\partial W}  & \frac{\partial^2 f}{\partial A^2}
\end{bmatrix}
\end{gather}
The Hessian matrix will be approximated by finite difference method (quasi Newton method), thus there is no need for detail analytical expression. The Hessian matrix has 6 non-zero elements.
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Model Predictive Controller Design 2}
\subsection{Mathematical formulation}
MPC formulation similar to MPC\#1 shown above, but only use $W$ and $A$ as manipulated variables:
\begin{gather}\label{eq:MPC}
\dis{\min_{W,A}} f(W,A) = q_{r^2}\left(r^2_{set}-\left<r^2_f\right>\right)^2+q_{m^2}\left(m^2_{set}-\left<m^2_f\right>\right)^2
\end{gather}
%\operatornamewithlimits{exp}_{x \in \mathbb{R}^n} J(t) = q_{r^2}\left(r^2(t_f)-r^2_{set}\right)^2+q_{m^2}\left(m^2(t_f)-m^2_{set}\right)^2\\
where
\begin{align}
\left<r^2_f\right> = \frac{1}{L^2}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{\infty}\sum_{p=1}^{4}\left<z^2_{p,n_x,n_y}(t_f)\right>,\quad
\left<m^2_f\right> = \sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}\sum_{p=1}^{4}\left(K_{p,n_x,n_y}\left<z^2_{p,n_x,n_y}(t_f)\right>\right)
\end{align}
\begin{align}
\left<z^2_{p,n_x,n_y}(t_f)\right> &=\var(z_{p,n_x,n_y}(t_f))+\left<z_{p,n_x,n_y}(t_f)\right>^2\\
\left<z_{p,n_x,n_y}(t_f)\right>   &= e^{\lambda_{n_x,n_y}(t_f-t)}\left<z_{p,n_x,n_y}(t)\right>+\frac{w_p}{\lambda_{n_x,n_y}}(e^{\lambda_{n_x,n_y}(t_f-t)}-1)\\
\var(z_{p,n_x,n_y}(t_f))     &= e^{2\lambda_{n_x,n_y}(t_f-t)}\var(z_{p,n_x,n_y}(t))+\sigma^2\frac{e^{2\lambda_{n_x,n_y}(t_f-t)}-1}{2\lambda_{n_x,n_y}}\\
\lambda_{n_x,n_y}            &= -\frac{4c_2\pi^2}{L^2}(n_x^2+n_y^2), \quad n_x^2+n_y^2\neq0
\end{align}
\begin{align}
%c(T,W)       &= W\left(1-\frac{k_w}{W^{a_w}}e^{k_BT/E_w}\right)\\
c_2(W)      &= p_{c1}W^4+p_{c2}W^3+p_{c3}W^2+p_{c4}W+p_{c5}\\
\sigma^2(W) &= p_{s1}W+p_{s2}
\end{align}
Subject to:
\begin{gather}
W_{min}\leq W \leq W_{max}, \quad |W(t)-W(t-dt)|\leq \Delta W_{max},\label{eq:MPC_LastConstraint}
\end{gather}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Implementation using the IPOPT optimizer}
In order to implement in IPOPT, the optimization problem in eq.~\eqref{eq:MPC}--\eqref{eq:MPC_LastConstraint} is rewritten as
\begin{itemize}
\item{
Cost function
\begin{gather}\label{eq:MPC}
f(W,A) = q_{r^2}\left(r^2_{set}-\left<r^2_f\right>\right)^2+q_{m^2}\left(m^2_{set}-\left<m^2_f\right>\right)^2
\end{gather}
%\operatornamewithlimits{exp}_{x \in \mathbb{R}^n} J(t) = q_{r^2}\left(r^2(t_f)-r^2_{set}\right)^2+q_{m^2}\left(m^2(t_f)-m^2_{set}\right)^2\\
where
\begin{align}
\left<r^2_f\right> = \frac{1}{L^2}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{\infty}\sum_{p=1}^{4}\left<z^2_{p,n_x,n_y}(t_f)\right>,\quad
\left<m^2_f\right> = \sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}\sum_{p=1}^{4}\left(K_{p,n_x,n_y}\left<z^2_{p,n_x,n_y}(t_f)\right>\right)
\end{align}
\begin{align}
\left<z^2_{p,n_x,n_y}(t_f)\right> &=\var(z_{p,n_x,n_y}(t_f))+\left<z_{p,n_x,n_y}(t_f)\right>^2\\
\left<z_{p,n_x,n_y}(t_f)\right>   &= e^{\lambda_{n_x,n_y}(t_f-t)}\left<z_{p,n_x,n_y}(t)\right>+\frac{w_p}{\lambda_{n_x,n_y}}(e^{\lambda_{n_x,n_y}(t_f-t)}-1)\\
\var(z_{p,n_x,n_y}(t_f))     &= e^{2\lambda_{n_x,n_y}(t_f-t)}\var(z_{p,n_x,n_y}(t))+\sigma^2\frac{e^{2\lambda_{n_x,n_y}(t_f-t)}-1}{2\lambda_{n_x,n_y}}\\
\lambda_{n_x,n_y}            &= -\frac{4c_2\pi^2}{L^2}(n_x^2+n_y^2), \quad n_x^2+n_y^2\neq0
\end{align}
\begin{align}
%c(T,W)        &= W\left(1-\frac{k_w}{W^{a_w}}e^{k_BT/E_w}\right)\\
c_2(W)      &= p_{c1}W^4+p_{c2}W^3+p_{c3}W^2+p_{c4}W+p_{c5}\\
\sigma^2(W) &= p_{s1}W+p_{s2}
\end{align}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Gradient of cost function
\begin{gather}
\nabla f = \left[\frac{\partial f}{\partial W},\frac{\partial f}{\partial A}\right]
\end{gather}
\begin{gather}
\begin{split}
\frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial \lambda_{n_x,n_y}} =& 2e^{2\lambda_{n_x,n_y}(t_f-t)}\var(z_{p,n_x,n_y}(t))(t_f-t)\\
                                                                      &+\sigma^2\frac{4\lambda_{n_x,n_y}e^{2\lambda_{n_x,n_y}(t_f-t)}(t_f-t)-2\left(e^{2\lambda_{n_x,n_y}(t_f-t)}-1\right)}{4\lambda_{n_x,n_y}^2}
\end{split}
\end{gather}
\begin{gather}
\frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial \sigma^2} = \frac{e^{2\lambda_{n_x,n_y}(t_f-t)}-1}{2\lambda_{n_x,n_y}}
\end{gather}
\begin{gather}
\begin{split}
\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial \lambda_{n_x,n_y}} &=e^{\lambda_{n_x,n_y}(t_f-t)}\left<z_{p,n_x,n_y}(t)\right>(t_f-t)\frac{\partial \lambda_{n_x,n_y}}{\partial T}\\
&+w_p\frac{\lambda_{n_x,n_y}e^{\lambda_{n_x,n_y}(t_f-t)}(t_f-t)-(e^{\lambda_{n_x,n_y}(t_f-t)}-1)}{\lambda^2_{n_x,n_y}}
\end{split}\\
\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial \sigma^2}=0
\end{gather}
%%%%%%%%%%%%%%%%%%% \frac{\partial f}{\partial W} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
$\frac{\partial f}{\partial W}$ is given by
\begin{gather}
\frac{\partial f}{\partial W}= -2q_{r^2}\left(r^2_{set}-r^2_f\right)\frac{\partial r^2_f}{\partial W}-2q_{m^2}\left(m^2_{set}-m^2_f\right)\frac{\partial m^2_f}{\partial W}
\end{gather}
\begin{gather}
\frac{\partial r^2_f}{\partial W} = \frac{1}{L^2}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{\infty}\sum_{p=1}^{4}\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial W}\\
\frac{\partial m^2_f}{\partial W} = \sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{\infty}\sum_{p=1}^{4}\left(K_{p,n_x,n_y}\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial W}\right)
\end{gather}
\begin{gather}
\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial W}=\frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial W}+2\left<z_{p,n_x,n_y}\right>\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial W}
\end{gather}
\begin{gather}
\frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial W} = \frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial \lambda_{n_x,n_y}}\frac{\partial \lambda_{n_x,n_y}}{\partial W}+\frac{\partial \var(z_{p,n_x,n_y}(t_f))}{\partial \sigma^2}\frac{\partial \sigma^2}{\partial W}
\end{gather}
\begin{gather}
\frac{\partial \lambda_{n_x,n_y}}{\partial W}= -\frac{4\pi^2}{L^2}(n_x^2+n_y^2)\frac{\partial c_2}{\partial W}\\
\frac{\partial c_2}{\partial W} = 4p_{c1}W^3+3p_{c2}W^2+2p_{c3}W+p_{c4}
\end{gather}
\begin{gather}
\frac{\partial \sigma^2}{\partial W} = p_{s1}
\end{gather}
\begin{gather}
\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial W}=\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial \lambda_{n_x,n_y}}\frac{\partial \lambda_{n_x,n_y}}{\partial W}+\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial w_{p,m,n}}\frac{\partial w_{p,m,n}}{\partial W}\\
\frac{\partial w_{p,m,n}}{\partial W}=
\begin{cases}
L    & p=2,m,n=0\\
0    & \text{otherwise}
\end{cases}
\end{gather}
%%%%%%%%%%%%%%%%%%% \frac{\partial f}{\partial A} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
$\frac{\partial f}{\partial A}$ is given by
\begin{gather}
\frac{\partial f}{\partial A}= -2q_{r^2}\left(r^2_{set}-r^2_f\right)\frac{\partial r^2_f}{\partial A}-2q_{m^2}\left(m^2_{set}-m^2_f\right)\frac{\partial m^2_f}{\partial A}
\end{gather}
\begin{gather}
\frac{\partial r^2_f}{\partial A} = \frac{1}{L^2}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{\infty}\sum_{p=1}^{4}\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial A}\\
\frac{\partial m^2_f}{\partial A} = \sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq 0}}^{\infty}\sum_{p=1}^{4}\left(K_{p,n_x,n_y}\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial A}\right)
\end{gather}
\begin{gather}
\frac{\partial \left<z^2_{p,n_x,n_y}(t_f)\right>}{\partial A}=2\left<z_{p,n_x,n_y}\right>\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial A}
\end{gather}
\begin{gather}
\begin{split}
\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial A}=\frac{\partial \left<z_{p,n_x,n_y}(t_f)\right>}{\partial w_p}\frac{\partial w_p}{\partial A} = \frac{1}{\lambda_{n_x,n_y}}\left(e^{\lambda_{n_x,n_y}(t_f-t)}-1\right)\frac{\partial w_{p,n_x,n_y}}{\partial A}
\end{split}\\
\frac{\partial w_{p,n_x,n_y}}{\partial A} =
\begin{cases}
0                                         & \text{if $p=1$}\\
\frac{L}{2k\pi}\left[1-\cos(2k\pi)\right] & \text{if $p=2, n_x,n_y=0$}\\
0                                         & \text{if $p=2, n_x,n_y \neq 0$}\\
0                                         & \text{if $p=2, n_x=k, n_y  = 0$}\\
\frac{\sqrt{2}Lk}{2\pi(n_x^2-k^2)}\left[\cos(2k\pi)-1\right] & \text{if $p=2, n_x\neq k, n_y = 0$}\\
0                                         & \text{if $p=3, n_x,n_y = 0$}\\
0                                         & \text{if $p=3, n_x,n_y \neq 0$}\\
\frac{\sqrt{2}L}{2}                       & \text{if $p=3, n_x=k, n_y =0$}\\
\frac{\sqrt{2}Ln_x}{2\pi(k^2-n_x^2)\sin(2k\pi)} & \text{if $p=3, n_x\neq k, n_y=0$}\\
0                                         & \text{if $p=4$}
\end{cases}
\end{gather}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Inequality constraints
\begin{gather}
0\leq g(x) = W-A \leq \infty
\end{gather}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Jacobian matrix of constraints
\begin{gather}
\nabla g = \begin{bmatrix} 1 & -1 \end{bmatrix}
\end{gather}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Variable constraints
\begin{gather}
\max\{W_{min},W_{old}-\Delta W_{max}\}\leq W \leq \min\{W_{max},W_{old}+\Delta W_{max}\}\\
0 \leq A \leq \infty
\end{gather}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Hessian matrix of the Lagrangian
\begin{gather}
\sigma_f\nabla^2 f+\lambda_1\nabla^2 g=
\sigma_f
\begin{bmatrix}
\frac{\partial^2 f}{\partial W^2}          & \\
\frac{\partial^2 f}{\partial A\partial T}  & \frac{\partial^2 f}{\partial A^2}
\end{bmatrix}
\end{gather}
The Hessian matrix will be approximated by finite difference method (quasi Newton method), thus there is no need for detail analytical expression. The Hessian matrix has 3 non-zero elements.
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{itemize}
\end{document}
